^  CO 

e-  CO; 


^  r 

h:  '"~:J* 

^rTl  !&*.»*.  *** 


Q 

«£ 


*BC  OO 
O  <3C 


^3-  3-3 


Office  of  Naval  Research 
Department  of  the  Navy 

Contract  Nonr-220(28) 


403  347 


THE  COLLAPSE  OF  A  SPHERICAL  CAVITY 
IN  A  COMPRESSIBLE  LIQUID 


Robert  Hickling  and  Milton  S.  Plesset 


Division  of  Engineering  and  Applied  Science 
CALIFORNIA  INSTITUTE  OF  TECHNOLOGY 
Pasadena,  California 


Report  No.  85-24 


March,  1963 


Office  of  Naval  Research 
Department  of  the  Navy 
Contract  Nonr- 220(28) 


THE  COLLAPSE  OF  A  SPHERICAL  CAVITY 
IN  A  COMPRESSIBLE  LIQUID 

Robert  Hickling  and  Milton  S.  Plesset 


Reproduction  in  whole  or  in  part  is  permitted  for  any  purpose  of 
the  United  States  Government 


Division  of  Engineering  and  Applied  Science 
California  Institute  of  Technology 
Pasadena,  California 


Report  No.  85-24 


March,  1963 


THE  COLLAPSE  OF  A  SPHERICAL  CAVITY 
IN  A  COMPRESSIBLE  LIQUID 

by 

Robert  Hickling  and  Milton  S.  Plesset 
Abstract 

This  paper  presents  numerical  solutions  for  the  flow  in  the 
vicinity  of  a  collapsing  spherical  bubble  in  water.  The  bubble  is  assumed 
to  contain  a  small  amount  of  gas  and  the  solutions  are  taken  beyond  the 
point  where  the  bubble  reaches  its  minimum  radius  up  to  the  stage  where 
a  pressure  wave  forms  and  propagates  outwards  into  the  liquid.  The 
motion  up  to  the  point  where  the  minimum  radius  is  attained,  is  found  by 
solving  the  equations  of  motion  both  in  the  Lagrangian  and  in  the  char¬ 
acteristic  forms.  These  are  in  good  agreement  with  each  other  and  also 
with  the  approximate  theory  of  Gilmore  which  is  demonstrated  to  be 
accurate  over  a  wide  range  of  Mach  number.  The  liquid  flow  after  the 
minimum  radius  has  been  attained  is  determined  from  a  solution  of  the 
Lagrangian  equations.  It  is  shown  that  an  acoustic  approximation  is 
quite  valid  for  fairly  high  pressures  and  this  fact  is  used  to  determine 
the  peak  intensity  of  the  pressure  wave  at  a  distance  from  the  center  of 
collapse.  It  is  estimated  in  the  case  of  typical  cavitation  bubbles  that 
such  intensities  are  sufficient  to  cause  cavitation  damage. 
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1 .  Introduction 

The  effect  of  compressibility  in  the  collapse  of  a  cavity  in  a  liquid 
has  been  studied  for  some  years.  Rayleigh^  ^  was  the  first  to  point  out 
that  this  was  necessary  when  he  showed  that  pressures  in  the  liquid 
adjacent  to  the  cavity  wall  could  reach  very  high  values.  However, a  real 
interest  in  the  problem  originated  some  time  later,  out  of  work  on  under¬ 
water  explosions ^  .  Analytical  theories  were  developed  by  Herring^  ^ 

[4] 

and  Trilling1,  J  which  were  first  order  approximations  based  oi* 

Rayleigh's  solution^1  ^  for  an  incompressible  liquid.  To  obtain  higher 

order  approximations  the  so-called  Kirkwood -Bethe  hypothesis^  ^ 

was  postulated.  The  Kirkwood-Bethe  hypothesis  forms  the  basis  of  a 

solution  by  Gilmore^  ^  which  has  been  used  in  a  modified  form  by 
r7i  r  3 1 

Mellen1  1  and  by  Flynn1  .  Gilmore's  solution  has  been  found  to  be 
surprisingly  accurate  when  compared  to  exact  solutions  obtained  numerical¬ 
ly,  and  this  behavior  is  again  confirmed  by  our  results. 

With  the  advent  of  high-speed  computers  numerical  solutions  of 
the  exact  equations  of  compressible  flow  became  possible,  and  several 

integrations  of  the  characteristic  equations  have  been  performed. 

[91 

Gilmore  1  obtained  some  results  which  he  used  to  compare  with  his 
approximate  theory,  while  Hunter^**^  performed  some  numerical  work 
which  acted  as  a  basis  for  a  similarity  solution  valid  in  the  immediate 
neighborhood  of  the  collapse  point  of  an  empty  cavity.  Brand^1^  did 
some  further  computations.  So  far  however,  there  has  been  no  numerical 
solution  which  describes  the  rebound  of  the  cavity  and  the  subsequent 
formation  of  a  shock  wave  traveling  outwards  into  the  liquid.  The  work 
of  this  paper  is  concerned  with  such  a  solution. 
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Classically  the  problem  to  be  solved  is  that  of  the  empty  cavity. 

This  is  usually  quite  a  good  model  since  the  small  amounts  of  gas,or 
vapor,  which  occur  in  a  typical  cavitation  bubble,  have  little  effect  on 
the  motion  of  the  interface  until  the  final  stage  of  the  collapse.  However, 
the  contents  of  a  cavity,  even  though  they  may  be  quite  rarified  initially, 
will  have  an  important  effect  on  the  final  pressure,  since  the  further  a 
collapse  can  proceed  before  being  halted  the  greater  the  pressure.  The 
final  stage  of  the  collapse  will  be  affected  by  the  behavior  of  the  gas  in¬ 
side  the  cavity  as  well  as  by  the  quantity.  Thus  an  isothermal  gas  can  be 
expected  to  produce  a  more  violent  collapse  t'  ~>n  an  equal  amount  of 
adiabatic  gas.  It  is  not  our  purpose  in  this  paper  to  consider 
the  physical  behavior  at  the  final  collapse  in  detail.  Such  a  task 
would  in  fact  be  quite  difficult.  The  main  interest  here  is  in  the  effect 
of  the  liquid  compressibility.  Thus  the  behavior  of  the  contents  of  the 
cavity  will  merely  be  approximated  by  a  simple  gas -like  model.  The  cavity 
will  be  .assumed  to  contain  a  uniform  gas  whose  pressure  varies  accord¬ 
ing  to  the  law  p~  ,  where  p  is  the  density.  The  index  y  can  be 
varied  so  as  to  simulate  the  gross  behavior  of  the  gas.  Thus  y  =  1 
implies  that  the  gas  is  isothermal  while  y  >  I  implies  that  some  of  the 
heat  of  compression  is  being  retained  in  the  gas.  This  model  has  the 
advantage  of  avoiding  the  singularities  which  would  occur  at  the  final 
collapse  point  if  the  cavity  were  empty.  The  sudden  compression  of  a 
gas  or  vapor  inside  a  cavity  can  be  expected  to  raise  its  temperature, 
and  due  to  conduction  and  condensation  there  will  be  a  consequent  rise 

in  temperature  in  the  liquid  at  the  interface.  These  effects  will  be 

F 1 2 1 

confined  to  a  very  narrow  region1  1 ,  and  they  will  therefore  be  neglected. 
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A  sudden  compression  of  a  liquid  does  not  cause  a  significant  increase  in 
temperature^  *  ^ ,  so  that  it  will  be  assumed  to  behave  according  to  the 
well-known  law  pn«(p+B)  where  B  and  n  are  constants.  The  results 
presented  here  are  for  water  with  B  equal  to  3,000  atmospheres  and  the 
index  n  equal  to  7.  The  validity  of  these  values  has  been  discussed 
elsewhere^^'  It  is  also  assumed  that  the  liquid  is  inviscid  and 

that  the  motion  remains  spherically  symmetric  at  all  times.  This 
description  of  a  collapsing  cavity  is  of  course  not  completely  accurate  in 
detail,  but  it  is  considered  to  be  sufficiently  satisfactory  for  a  study  of 
the  effects  of  liquid  compression. 

Numerical  solutions  of  the  equations  of  compressible  flow  for  a 
collapsing  cavity  in  water  are  presented  here  for  both  the  characteristic 
and  the  Lagrangian  forms  of  the  equations.  A  comparison  between  the 
two  provides  a  check  on  the  numerical  accuracy.  The  Lagrangian 
solution  was  carried  beyond  the  final  collapse  point  into  the  region  where 
the  liquid  rebounds  and  generates  a  compression  wave  which  travels  out¬ 
wards  from  the  collapse  center.  The  solution  was  terminated  when  the 
compression  wave  had  steepened  into  a  shock  front.  It  could  probably 
have  been  continued  beyond  this  stage  by  using  techniques  such  as  those 
proposed  by  von  Neumann  and  Richtmeyer^*^  and  Lax^*^,  where  an 
artificial  viscosity  is  introduced  into  the  equations  of  motion  to  produce 
a  continuous  variation  through  the  shock.  However,  enough  information 
had  already  been  obtained  to  estimate  the  order  of  magnitude  of  peak  in¬ 
tensities  at  a  distance  from  the  collapse  center  so  that  no  attempt  was 
made  to  undertake  this  additional  calculation. 

One  of  the  principal  purposes  of  this  analysis  was  to  determine 
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whether  shock#  emanating  from  the  collapse  of  cavitation  bubbles  could 
provide  a  mechanism  for  cavitation  damage.  The  results  given  here  show 
that  such  shocks  can  be  strong  enough  to  cause  damage  to  metals  and 
other  solids  in  the  vicinity  of  the  bubble. 

2.  Formulation  of  the  Problem 

A  spherical  cavity  containing  a  uniform  gas  is  assumed  to  expand 
or  contract  in  an  infinite  volume  of  liquid.  For  cavitation  bubbles  the 
effect  of  gravity  is  generally  quite  small,  so  that  it  will  be  neglected 
along  with  other  asymmetric  effects.  Since  the  motion  is  spherically 
symmetric,  it  will  be  irrotational  i.e. 

curl  u  =  0  (1) 

where  u  is  the  velocity  vector.  If  the  center  of  the  bubble  is  chosen  as 
the  center  of  a  system  of  spherical  polar  coordinates,  then  u  will  have 
only  a  single  component  u  lying  in  the  radial  direction.  The  equations 
of  motion  expressing  the  conservation  of  mass  and  momentum  for  a 
spherically  symmetric  system  are 

»-*[£♦£]  •  <*> 


where  p  and  p  are  the  pressure  and  density  in  the  liquid,  and  the 


operator 


D  8 
TH  ”  "5T 


+  u 
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■Sr 


(4) 
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is  the  derivative  with  respect  to  time  following  the  motion  of  the,  liquid. 
It  is  assumed  that  the  liquid  is  isentropic^  ^  and  that  the  density  and 
pressure  are  related  by  an  equation  of  state  of  the  form 


where  p^.  and  p^  are  the  pressure  and  density  in  the  liquid  at  infinity. 
For  water  the  constant  B  is  given  a  value  of  3,000  atmospheres  while 
the  index  n  has  the  value  7.  The  use  of  this  equation  has  been  justified 
by  several  authors  ^  ^  ^  ^  and  is  based  on  the  fact  that  entropy  changes 
are  small  even  when  very  large  pressure  jumps  are  present.  The  upper 
limit  on  the  accuracy  of  the  formula  appears  to  be  for  pressures  of  about 

5 

10  atmospheres.  The  velocity  of  sound  c  in  the  liquid  is  defined  by 

c*  =  =  l!iE±5l  .  (6) 

Op  p 

The  velocity  of  sound  at  infinity  is  therefore  given  by  cM=  [  n(pQo+B)/p(jc]^  . 

The  equations  of  motion  given  above  are  in  the  Eulerian  form 
which  describes  what  is  happening  at  a  particular  point  in  space  and  how 
it  varies  with  time.  However,  it  is  preferable  here  to  use  the  Lagrangian 
form  where  the  properties  of  the  fluid  are  obtained  by  following  the  par¬ 
ticle  motion.  To  do  this,  it  is  necessary  to  ascribe  to  each  particle  of  fluid 
a  value  of  a  certain  parameter,  y,  which  is  defined  by  the  expression 

r  r(y,t) 

y  =  \  pr‘dr  . 

Jr(o,t) 


Hence 
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P  r 


dr 

"®y 


r  1 


(7) 


With  this  relation,  Eq. (4)  becomes 

but,  since  Eq.(7)is  a  solution  of  this  equation,  the  two  conditions  are 
equivalent. 

On  the  walls  of  the  cavity  the  pressure  is  given  by  the  relation 

P  =  P0(R0/R)3  Y  (9) 

where  pQ  is  the  initial  pressure  inside  the  cavity  prior  to  collapse  and 
Rq  is  the  initial  radius  of  the  cavity.  Capital  letters  are  used  to  denote 
the  values  of  variables  at  the  cavity  wall.  Thus  R  is  the  radius  of  the 
cavity  and  U  its  velocity.  At  infinity  the  liquid  is  at  rest  and  the  pres¬ 
sure  and  density  have  the  values  p^  and  p^.  All  that  now  remains  be¬ 
fore  solving  the  equations  is  to  specify  the  pressure  and  velocity  in  the 
liquid  at  some  initial  instant.  For  a  cavity  in  an  incompressible  liquid, 
it  is  usual  to  suppose  that  the  liquid  is  initially  at  rest  and  at  a  uniform 
pressure  p%  and  that  the  collapse  is  generated  by  a  pressure  discon¬ 
tinuity  at  the  cavity  wall  i.e.  pw>  PQ.  The  instant  this  pressure  dis¬ 
continuity  is  allowed  to  take  effect,  it  is  felt  throughout  the  entire  volume 
of  the  liquid  because  of  the  incompressibility.  The  cavity  wall  then  starts 
to  accelerate  inward  from  rest.  For  a  compressible  liquid  there  would 
be  an  initial  jump  in  velocity  during  an  infinitesimal  period  of  time. 
Gilmore  has  shown  that  this  instantaneous  increase  in  velocity  is 
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given  by  the  relation 

uo  '  isrir  <«.*'«>  <10> 

where  c  -  [n(p  +B)/p  ]  .  The  maximum  value  of  D  ueed  here  ie 

O  O  O  0 G 

10  atmospheres,  so  that  for  water,  such  a  jump  in  velocity  will  always 
be  quite  small  compared  to  sonic  velocities.  In  fact,  during  the  early 
stages  of  the  collapse,  the  solution  will  be  indistinguishable  from  that  for 
the  incompressible  liquid.  Because  of  their  simplicity  the  same  initial 
conditions  will  be  used  in  the  present  analysis  and,  even  though  it  is  not 
really  necessary,  the  effect  of  the  initial  jump  in  velocity  will  be  in¬ 
cluded.  The  possibility  of  the  formation  of  an  initial  compression  wave 
in  the  gas  or  vapor  inside  the  cavity  is  neglected.  The  initial  conditions 
are  somewhat  artificial  and  there  is  not  much  to  be  gained  here  in  follow¬ 
ing  up  such  an  implication. 

A  method  widely  used  in  the  solution  of  the  equations  of  compres¬ 
sible  flow  is  the  method  of  characteristics.  The  characteristic  equations 
for  spherically  symmetric  flow  are 


9r 

"35 

=  (u+c) 

9t 

"55  * 

(ll> 

dr 

■sp 

=  (u-c) 

9t 

"3$  ’ 

(12) 

0U  , 

c  9p 

+  2cu 

dt  n 

S5  + 

p  "5a 

r 

"55  -°  ' 

8u 

c  80 

2cu 

8t  n 

P  "5p 

r 

"S?  =  0  ‘ 

The  last  two  equations  are  not  very  suitable  for  finite  difference  work 
since  p  varies  quite  slowly  in  the  liquid.  Since  6p=6p/c*,  the 
alternative  forms 


8u  ,  1  8p  ,  2cu  8t  n 
"3a  +  pc  Ta  +  T“  "85  "  0 

(13) 

8u  ^  1  8p  _  2cu  8t  0 

(14) 

5P  pc  "3p  r  3p 

i 
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are  used.  Equations  (1 1)  and  (1 2)  define  the  system  of  characteristic  lines 
a  and  0  where  a  is  the  outward-going  characteristic  and  p  the  in¬ 
ward  going.  The  bubble  wall  motion  is  defined  by  the  variable  I  and 
is  determined  from  the  expressions 


and 


8R 

W 


(15) 


(16) 


Related  to  the  method  of  characteristics  is  the  approximate  theory 
of  Gilmore^]  which  is  based  on  the  so-called  Kirkwood- Be  the  assump- 

r  y  1  2 

tion*-  .  This  states  that  the  quantity  r  £h(p)  +  j  is  a  constant  along 
an  outward  going  characteristic,  where 


h(p)  =  J  =  (C*'  C*>  U7) 

Poc 

is  the  enthalpy  difference  between  the  liquid  at  pressure  p  and  at  pres¬ 
sure  p^.  under  isentropic  conditions.  This  assumption  can  be  express¬ 
ed  in  the  form  of  an  equation  containing  the  particle  derivative  as  defined 
in  Eq.  (4).  Since  the  cavity  wall  moves  with  the  liquid,  and  since  this 
motion  can  be  expressed  purely  as  a  function  of  time,  such  particle 
derivatives  can  be  changed  into  ordinary  derivatives  with  respect  to  time 
and  the  equation  will  become 


R 


dU 

cTT 


r,  u  i 

,  3  TT*  r 

in  „  r 

1  X  U1 

X  R 

dH  r 

1  U 1 

L1  '  z\ 

+ 1  u  1 

3cJ  "H  1 

+  c 

HF  L 

1  ’  -cl 

(18) 


which  governs  the  motion  of  the  cavity  wall .  Given  the  initial  values  on  the 
boundary,  a  second  ordinary  differential  equation  gives  the  liquid  con¬ 
ditions  along  an  outward-going  characteristic  line: 
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du  _  R(H+U2  /2)(u+c) 
^  r2(c-u) 


2c2u 
7{cmI)  • 


(19) 


This  is  solved  in  conjunction  with  the  Kirkwood-Bethe  assumption  and 
Eqs.  (ll)and(17). 

For  most  problems  of  interest  in  bubble  collapse  it  can  be 
assumed  that  |  H |  <  <  C2.  With  water  for  example,  this  inequality 
corresponds  to  a  pressure  difference  |  Pq^P^I  <  <  20,000  atmos.  Using 
this  assumption,  Gilmore was  able  to  derive  a  simple  relation  be¬ 
tween  the  bubble  wall  velocity  and  the  radius,  for  the  case  of  an  empty 
cavity.  This  relation  is 


3P«U*  1 


(20) 


It  is  seen  that  as  R  tends  to  zero,  U  varies  as  R 


-i 


in  contrast  to 


Til  -  yz 

incompressible  theory1  J  which  gives  U  varying  as  R 

The  variables  used  in  the  above  equations  were  non-dimensional- 

ized  with  respect  to  the  initial  radius  Rq,  the  density  p^.,  and  velocity 

of  sound  c  in  the  following  way: 
oo 


p  =  p/p*  ;  p'  =  p/p«c«  J 

u'  =  u/c 

00 

t'  =  c  t/R  ;  h'  =  h/c  2  ; 

00  O  00 

r’  = 

It  will  be  supposed  that  all  the  variables  and  constants  used  above  have 
already  been  non-dimensionalized  in  this  way,  and  the  primes  suppres¬ 
sed.  The  equations  of  motion  will  remain  unaltered  by  this  change,  but 
the  liquid  equation  of  state  Eq.(5)becomes 
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n(p  +  B)  =  pn  ,  (21) 

and  the  pressure  at  the  bubble  wall  given  by  Eq.  (9)becomes 

P=PqRjY  .  (22) 

Equations  (10)and  (1 7)are  modified  by  setting  ca=  1. 

It  will  be  seen  from  this  non-dimensionalisation  that  the  solutions 
are  actually  independent  of  the  initial  radius  Rq,  i.e. ,  that  the  same 
pressures  and  velocities  are  obtained  regardless  of  the  scale.  (The 
elapsed  time  is  of  course  proportional  to  Rq).  This  is  due  to  the  fact 
that  compressible  flow  equations  contain  only  first  order  derivatives. 

If  effects  such  as  heat- conduction ,  viscosity  and  surface  tension  were 
included,  the  solutions  would  become  dependent  on  RQ,  the  dependence 
becoming  stronger  with  a  decrease  in  RQ. 

Solutions  to  the  problem  of  bubble  collapse  were  obtained  by 
means  of  a  high-speed  computer,  using  the  three  methods  given  above. 
First  of  all  Gilmore's  method  was  used  to  establish  an  approximate 
solution  for  the  motion  of  the  cavity  wall.  The  initial  conditions  given 
above  were  applied,  together  with  the  initial  value  of  the  velocity  given 
by  Eq.  (10).  The  validity  of  this  solution  was  checked  in  relation  to  simple 
incompressible  flow  theory^ 1  ^  for  the  early  stages  of  the  motion.  Some 
exploratory  calculations  were  then  carried  out  using  the  method  of 
characteristics  in  the  region  where  the  Mach  number  U/C  at  the  cavity 
wall  was  of  the  order  of  0. 1.  These  results  were  compared,  both  on  the 
wall  and  in  the  interior  of  the  liquid,  with  results  obtained  by  Gilmore's 
method,  and  it  was  established  that  the  method  was  quite  accurate  in 
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this  region.  This  result  was  of  course  to  be  expected  since  it  can  be 
shown  ^  that  Gilmore's  theory  is  accurate  to  terms  at  least  of  the  order 
of  (u/c)a.  Thereafter,  Gilmore's  method  was  used  to  provide  initial 
values  in  the  subsonic  (U/C  ~0.  1)  region  for  the  solution  of  the 
Lagrangian  and  the  characteristic  equations.  Such  a  procedure  was 
necessary  because  of  the  large  amount  of  computing  time  required  to 
perform  an  exact  solution  starting  from  the  initial  stages  of  the  motion. 

The  Lagrangian  and  the  characteristic  solutions  were  carried  to  the  final 
collapse  point  and,  since  they  are  both  supposed  to  be  exact,  they  pro-1 
vided  a  check  on  the  numerical  accuracy  up  to  this  stage.  Comparisons 
were  made  for  points  on  the  cavity  wall  and  in  its  vicinity,  and  the 
discrepancies  which  occurred  were  estimated  to  be  at  most  about  one 
per  cent.  The  Lagrangian  solution  was  then  carried  on  into  the  region  of 
rebound  up  to  the  point  where  a  shock  wave  formed.  Solutions  were  ob¬ 
tained  for  a  variety  of  conditions  such  as  might  occur  in  cavitation.  The 
values  used  for  the  ambient  pressure  in  the  liquid  were  1  atmosphere, 

and  10  atmospheres.  The  initial  pressure  pQ  in  the  gas  was  varied 
-1  -4 

from  10  atmospheres  to  10  atmospheres,  while  the  index  y  for  the 
gas  had  the  values  1  and  1.4.  The  problems  were  programmed  for  an 
IBM  7090  and  the  numerical  methods  are  described  in  the  appendices. 

3.  Calculated  Results 

The  first  set  of  calculations  were  concerned  with  the  motion  of 
the  bubble  wall  and  with  the  gross  effect  of  the  gas  content  of  the  bubble. 
Figures  1-3  show  the  variation  of  the  bubble  wall  Mach  number  with 
the  radius  under  a  variety  of  conditions.  The  amount  of  gas  is  determined 
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by  the  initial  pressure  pQ,  and  the  behavior  of  the  gas  is  varied  by 
changing  the  index  y.  Figure  1  shows  the  collapse  of  a  bubble  for 
different  initial  pressures  pQ  with  y  =  1.4,  acting  under  an  external 
pressure  p#  of  1  atmosphere.  It  is  seen  that  as  the  amount  of  gas 
diminishes,  the  motion  of  the  bubble  wall  becomes  more  and  more  rapid 
during  the  final  stages  of  the  collapse.  In  the  case  of  the  empty  cavity, 
the  velocity  increases  without  limit  as  the  bubble  grows  smaller.  The 
exact  solutions  are  compared  with  solutionai  derived  from  Gilmore's 
theory  and  from  the  theory  of  an  empty  cavity  in  an  incompressible 
liquid,  it  is  seen  that  the  compressibility  becomes  very  important 
during  the  final  stages  of  the  motion.  Figure  2  shows  the  result  of  in¬ 
creasing  the  external  pressure  p^  from  1  atmosphere  to  10  atmos¬ 
pheres.  A  comparison  with  Fig.  3  shows  that  Buch  an  increase  does  not 
affect  the  collapse  so  much  as  a  change  in  y  from  1 . 4  to  1 .  Thus  a 
gas  which  behaves  isotherm  ally  should  in  general  produce  a  more  violent 
collapse  than  an  adiabatic  gas  reacting  under  a  high  external  pressure. 

In  all  these  results,  it  iB  evident  that  the  predictions  of  Gilmore's  theory, 
based  on  the  Kirkwood-Bethe  assumption,  continue  to  be  surprisingly 
accurate. 

For  the  empty  cavity  under  1  atmosphere  external  pressure,  the 

0  785 

bubble  wall  velocity  tends  to  infinity  as  (RQ/R)  ’  .  The  value  of  the 

index  was  found  to  be  the  same  when  pw  was  10  atmospheres.  This 
result  is  in  good  agreement  with  a  similar  estimate  by  Hunter^*^.  By 
comparison,  Gilmore's  theory  predicts  a  rate  varying  as  (R0/R)0,5, 
while  the  incompressible  theory  gives  (RQ/R)*  ’  *. 

The  remaining  results  are  concerned  with  the  flow  in  the  liquid 
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around  the  collapsing  bubble  and  were  obtained  for  the  two  cases, 
pQ  =  10  and  10  atmospheres,  with  y  =  1.4  and  an  external  pres¬ 
sure  of  1  atmosphere.  In  both  these  cases  the  liquid  rebounds  and 
forms  a  compression  wave  which  moves  outwards  and  steepens  into  a 
shock  front.  The  occurrence  of  the  shock  front  causes  the  numerical 
solution  to  become  unstable  so  that  the  results  are  only  presented  up  to 
this  point.  Methods  are  available^ ^ for  overcoming  this  instabil¬ 
ity,  and  proceeding  with  the  solution.  But  this  continuation  would  have 
involved  an  additional  program  of  computation  and,  since  it  is  not  really 
necessary  for  the  purposes  of  this  paper,  it  was  not  pursued.  Figure  4  - 
7  show  the  distributions  of  Mach  number  and  pressure  in  the  liquid. 

These  are  given  for  successive  instants  in  time  which  is  expressed  in 
terms  of  t,  the  time  required  for  the  bubble  to  collapse  from  the  initial 

radius  R  to  the  final  minimum  radius.  The  formula  used  to  determine 
o 

4 

the  time  scale  is  (t  -  t)10  / T  where  t  is  the  time  elapsed  from  the 
start  of  the  motion.  The  collapse  time  T  was  determined  accurately 
from  the  numerical  solutions .  It  can  also  be  estimated  by  use  of  the 
incompressible  theory  for  an  empty  cavity  in  an  incompressible  liquid. 
The  incompressible  relation  gives 

t-0.91  R0(p„/pJ4  • 

Estimates  from  this  agree  to  within  less  than  1%  with  values  obtained  in 
the  calculated  examples. 

Figures  5(b)  and  7(b)  show  the  pressure  wave  forming  and  travel¬ 
ing  outwards  into  the  liquid.  Because  of  the  compressiblity,  the  change 
in  the  direction  of  motion  of  the  interface  is  communicated  to  the  liquid 
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nly  by  the  passage  of  the  pressure  wave,  and  hence  a  reversal  of  flow 

occurs  through  it,  as  shown  in  Figs.  4  and  6.  The  shock  wave  forms 

.4 

fairly  rapidly  when  the  initial  pressure  is  10  atmospheres,  because 

the  final  collapse  pressures  are  high.  In  the  case  where  the  initial  pres- 
_  3 

sure  is  10  atmospheres,  the  pressure  front  does  not  steepen  so  quickly. 
When  po  =  10” 2  atmospheres,  the  final  pressure  at  the  cavity  wall  is 

3 

seen,  from  Fig.  1,  to  be  about  10  atmospheres.  For  liquid  pressures 
of  this  magnitude,  the  compression  wave  will  steepen  only  very  slowly 
and  the  rate  with  which  it  steepens  is  reduced  as  it  moves  outwards,  be¬ 
cause  of  geometric  attenuation.  For  such  cases  the  compression  will 
behave  like  an  acoustic  wave,  which  does  not  alter  much  in  form  as  it 
propagates.  For  larger  amounts  of  gas  inside  the  bubble,  the  intensity 
of  the  compression  will  diminish  even  further  until  only  a  very  weak 
pulse  emanates  from  the  collapse. 

Figure  5(b)  shows  that  the  acoustic  approximation  is  reasonably 
valid  even  in  the  case  where  a  shock  front  develops.  The  pressure  front 
gradually  steepens,  but  the  wave  remains  of  approximately  the  same 
thickness  and  form  and  attenuates  as  1/r  as  it  propagates  outwards 
from  the  collapse  center.  The  last  stage  of  the  calculations  shows  that 
the  peak  pressure  in  the  wave  is  about  1,000  atmospheres  at  R./Rq  0.  3. 

Beyond  this  point  not  much  change  will  occur  due  to  dispersion,  and  the 

[21 

losses  due  to  entropy  changes  through  the  shock  will  be  negligible1  . 
Hence  the  pressure  pulse  will  propagate  outwards  like  an  acoustic  wave, 
and  at  R/Ro  ^  2  the  peak  intensity  will  be  of  the  order  of  two  hundred 
atmospheres.  For  smaller  amounts  of  gas,  the  peak  pressures  will  be 
larger  although  they  are  a  little  more  susceptible  to  the  effects  of 
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dispersion.  Figure  7(b),  for  example,  shows  that  the  peak  pressure  still 
has  an  approximate  1/r  dependence  giving  an  intensity  of  the  order  of 
1,000  atmospheres  when  R/Rq  ~  2.  A  much  stronger  effect,  however, 
would  result  if  the  gas  inside  the  bubble  were  to  lose  some  of  its  heat  of 
compression.  Figures  1  and  3  show  that  if  the  gas  is  isothermal,  it 
collapses  to  a  much  smaller  radius  for  the  same  initial  amount  of  gas 
when  compared  to  the  adiabatic  condition.  It  has  been  shown^^for 
values  of  Rq  <  1  mm,  that  heat  losses  from  the  gas  into  the  liquid  are 
quite  significant.  For  bubbles  of  this  size  therefore,  a  more  violent 
collapse  would  occur,  and  pressure  wave  intensities  of  several  thousand 
atmospheres  could  be  expected. 

4.  Application  of  Results  to  the  Theory  of  Cavitation  Damage 

A  well-known  consequence  of  the  formation  of  cavitation  bubbles 
in  a  liquid  is  the  damage  which  can  occur  to  adjacent  solid  surfaces. 

Such  damage  has  been  shown  to  be  largely  mechanical  in  origin,  al¬ 
though  the  precise  mechanism  is  as  yet  undecided.  It  has  been  suggested 
for  example  that  the  damage  is  caused  by  bubbles  adhering  to  the  solid 
and  collapsing  asymmetrically  in  such  a  way  that  jets  of  liquid  form  and 
strike  the  surface  at  high  speed.  There  is  some  indication^ 18 ^  that 
such  a  mechanism  might  operate.  The  results  of  the  present  paper 
show  that  damage  could  certainly  result  from  pressure  waves  originating 
from  bubbles  situated  at  a  short  distance  from  the  surface.  This  would 
avoid  the  requirement  that  the  destructive  bubbles  be  only  those  which 
are  adhering  to  the  surface. 

Compared  to  usual  kinds  of  material  damage,  cavitation  damage 

-2 

is  very  localized.  Cavitation  bubbles  are  usually  of  the  order  of  10  cm. 


17 


in  sise  i.e.  of  the  order  of  the  grain  else  in  a  typical  metal.  Hence  the 

attack  from  the  cavitation  bubbles  will  be  directed  towards  the  individual 

grains  rather  than  towards  a  large  group  of  them.  This  would  imply  that 

cavitation  stresses  only  need  to  be  of  the  order  of  die  yield  stresses  of 

the  individual  crystals.  In  addition  the  fact  that  the  damage  usually 

results  after  an  exposure  to  cavitation  of  a  certain  duration,  indicates 

the  existence  of  some  kind  of  fatigue  process.  The  cold  working  which 
[17] 

has  been  observed  during  the  early  stages  of  the  exposure  lends  sup¬ 
port  to  this  conclusion.  It  can  be  estimated  therefore  that  pressures  of 
the  order  of  two  or  three  thousand  atmospheres  should  be  sufficient  to 
cause  most  of  the  kinds  of  cavitation  damage  which  have  been  observed. 

It  is  possible  of  course  that  local  stress  concentrations  may  occur  due  to 
irregularities  in  surface  finish  and  grain  structure  and  hence  the  original 
pressures  need  not  be  as  high  as  this.  However.for  present  purposes  it 
is  not  really  necessary  to  consider  such  possibilities. 

The  results  discussed  in  the  previous  section  indicate  that  pres¬ 
sure  waves  of  the  right  order  of  intensity  can  occur  from  collapsing 
cavitation  bubbles.  Such  bubbles  usually  contain  only  very  small  amounts 
of  vapor  and  gas,  and  are  usually  less  than  1  mm.  in  size  prior  to  collapse. 
Hence  the  collapse  process  should  be  violent  enough  to  produce  pressure 
waves  of  the  right  intensity.  This  estimate  is  based  on  the  assumption 

of  a  spherically  symmetric  collapse.  Both  the  presence  of  the  boundary 

r  i9i 

and  possible  instability  effects 1  J  during  the  final  stages  of  the  motion 
cam  upset  this  symmetry  and  presumably  reduce  the  intensity  of  the 
resulting  pressure  waves.  It  is  believed  that  the  results  for 
the  spherical  motion  give  a  sufficiently  good  indication  of  the  order  of 
magnitude  of  the  pressures  which  can  occur  in  practice. 
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APPENDIX  I  -  THE  KIRKWOOD-BETHE  SOLUTION 

From  the  theory  given  in  Section  2,  the  governing  equation  in  non- 
dimensionalized  form  for  the  motion  of  the  cavity  wall  are: 


gjr  =  [H(C+U)  +  £  ^  (C-U)  -  i  U*(3C-U) ]/R(C-U) 


dR 

*3t 


=  U 


^  =  -3yPU/R[n(P+B)]1/n 

p  =  pc/R3v 

C  =  [^P+B)/"'1^2” 


d.l) 

(1.2) 

(1.3) 

(1.4) 

(1.5) 


together  with  the  initial  conditions  at  t  =  0, 


R  =  1  ;  P  =  Pc  ;  C  =  C(pQ)  ; 

U  =  tc(P0)-U  5  H  =  J^T)  tC*(po)-l]  .  (1.6) 

These  equations  were  solved  at  intervals  At  by  a  simple  step  by  step 
iterative  procedure  starting  from  the  initial  values  Eq.(1.6).  A  first 
estimate  of  U  at  t  =  At  was  obtained  by  substituting  the  initial  values 
into  Eq.(I.  1).  This  was  averaged  with  the  initial  value  of  U  and  used  in 
Eq.  (I.  2)to  obtain  a  first  estimate  of  R  at  t  =  At.  The  values  of  R  were 
then  averaged  in  the  same  way  as  the  values  of  U  and  used  in  Eqs.(1.4) 
and(I.  5)to  get  correspondingly  averaged  values  of  P  and  C.  The  aver¬ 
aged  values  obtained  so  far  were  then  used  in  Eq.  (I.  3)to  obtain  an 
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averaged  form  of  dH/dt  and  a  firet  estimate  of  H  at  t  *  At.  Thie  pro¬ 
cedure  was  repeated  step  by  step  for  each  interval  At.  In  approaching 
the  collapse  point,  the  variation  with  time  becomes  quite  rapid,  and  it  is 
necessary  to  keep  reducing  the  interval  At  in  order  to  preserve  accuracy. 
This  was  done  by  repeatedly  carrying  the  solution  up  to  the  point  where 
U  becomes  positive,  then  returning  to  a  suitable  earlier  point  and  repeat¬ 
ing  the  last  stages  of  the  calculation  for  a  new  interval  At/5,  where  At 
was  the  previous  interval. 

The  solution  for  the  motion  of  the  cavity  wall  provides  the  initial 
values  for  the  following  set  of  equations,  which  can  be  solved  along  an 


out-going  characteristic  in  the  interior  of  the  liquid, 

^  .  [l£±Hl  Y  -  2c*u]  /r(c-u)  (1.7) 

=  [n(p+B)]1/n  [4c*u*  -  (ctP)Z  y]  /r(c-u)  (1.8) 

Tfc  =<u+c)  (1.9) 

C  =  [n(p+B)](n"1)/2n  .  (1. 10) 

The  constant  Y  is  derived  from  the  initial  conditions  and  is  given  by 

Y  =  R(H+U*  /  2)  (1.11) 


where  R,  H,  U  are  values  on  the  cavity  wall.  This  set  of  equations  was 
solved  in  the  manner  already  described. 


20 


APPENDIX  II  -  THE  SOLUTION  OF  THE 
CHARACTERISTIC  EQUATIONS 


The  non-dimensionalized  equations  for  the  characteristic  solution 

are 


8r 
■5 a 

=  £ 

(II.  I) 

dr 

■8F 

"  <U-C)  % 

(II.  2) 

da 

ss 

+  -L  ®E  + 

pc  da 

2cu 

r 

o 

II 

(II.  3) 

da 

W 

.  _L  ®P 
pc  Up 

2cu 

r 

o 

II 

(II.  4) 

n 

P 

=  n(n+B) 

(II.  5) 

c  = 

[n(p+B)](n_1)/2n 

(II.  6) 

The  method  used  to  solve  this  system  is  the  one  used  by  Gilmore^  ^  and 
later  by  Brand^1^.  The  initial  points  lie  along  an  out-going  characteristic 
line  a,  with  the  analytic  continuation  being  provided  at  the  intersection 
with  the  boundary  by  the  conditions 


OR 


=  U 


8t 

w 


(H.  7) 
(II.  8) 


This  is  somewhat  shorter  than  the  method  used  by  Hunter  who  took 
the  initial  values  to  lie  along  a  line  t  =  const.  In  the  present  calcula¬ 
tions,  the  initial  values  were  provided  by  Gilmore's  method. 
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The  Eqs.(II.  1  -  II.  4)and(II.7)(  II.8)were  expressed  in  finite  dif¬ 
ference  form  and  applied  according  to  the  scheme  shown  in  Fig.  8.  The 
initial  line  of  points  along  an  out-going  characteristic  are  designated  by 
the  letter  i.  The  point  (i,j)  represents  the  jth  point  on  this  line  start¬ 
ing  from  the  point  on  the  boundary.  The  continuation  is  provided  at  the 
boundary  at  the  intersection  of  the  inward-going  characteristic  through 
(i,2)  with  the  boundary.  This  yields  the  new  point  (i  +  1 ,  1).  The  out¬ 
ward-going  characteristic  through  this  new  point  and  the  inward-going 
characteristic  through  the  point  (i,  3)  intersect  at  the  next  new  point 
(i  +  1,  2).  By  repeating  this  process,  one  obtains  all  the  new  points  on 
the  line  (i+1).  It  is  seen  that  the  total  number  of  these  points  is  dimin¬ 
ished  by  1  compared  with  the  line  i.  The  kind  of  averaging  technique 
described  in  the  previous  appendix  was  also  used  here  to  improve  the 
accuracy. 

The  fact  that  one  point  was  lost  at  each  step  of  the  solution  proved 
to  be  a  problem,  because  it  was  difficult  to  decide  how  many  points  would 
be  needed  to  carry  the  solution  to  the  final  stages  of  the  collapse.  This 
was  remedied  by  adding  one  point  at  each  step  using  Gilmore's  method. 
This  was  always  done  when  the  Mach  number  was  small  so  that  the  new 
point  could  be  considered  to  have  an  accuracy  corresponding  to  the  other 
points  on  the  characteristic  line.  As  the  solution  approached  the  last 
stages  of  the  collapse,  the  outward-going  characteristics  began  to  extend 
into  the  region  of  rebound  i.e.  into  the  region  t  >  t  where  t  is  the 
time  of  collapse.  The  main  criterion  for  this  was  the  time  when  the 
velocity  began  to  turn  positive  at  the  end  of  the  line  of  points.  Whenever 
this  occurred,  the  affected  points  were  discarded.  In  this  way  only  the 
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points  needed  to  give  a  solution  up  to  the  final  stages  of  collapse  were 
preserved. 

The  number  of  points  used  was  also  determined  by  the  motion  of 
the  bubble  wall.  When  the  collapse  becomes  more  rapid,  more  points 
are  needed  to  ensure  the  proper  accuracy.  As  a  check  the  first  and 
final  estimates  of  the  velocity  and  pressure  were  compared  at  each  step 
and,  if  these  differed  by  a  certain  percentage,  then  an  interpolation  auto¬ 
matically  occurred  which  doubled  the  total  number  of  points.  The  ad¬ 
ditional  points  were  situated  at  the  half  intervals. 

The  techniques  described  here  were  tested  for  a  variety  of 
examples,  and  the  accuracy  of  the  results  presented  is  considered  to  be 
quite  satisfactory.  The  calculations  were  performed  on  an  IBM  7090 
and  each  example  took  about  ten  minutes. 
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APPENDIX  m  -  SOLUTION  OF  THE  LAGRANGIAN  EQUATIONS 


are, 


The  governing  equations  for  the  Lagrangian  form  of  the  solution 


P* 


z 


dr 

W 


=  1 


(HI.  1 ) 


8u 

“St 


(HI.  2) 


with 


8r 

“5t 


=  u 


p"  =  n(p+B) 


(HI.  3) 

(HI.  4) 


cS[n(p+B)fn'1)/2n  . 


(in.  5) 


This  set  of  equations  is  sufficient  to  determine  the  variables.  Equation 
(III.  1) can  be  used  to  find  p,  whence  Eq. (HI. 4)yields  p.  However  this 
is  not  a  very  satisfactory  process,  since  p  varies  comparatively  slowly 
and  is  not  very  suitable  as  a  running  variable.  Hence  instead  of  Eq.(III.  1) 
the  following  relation  was  used, 

=  -  [n(p+B)]  [pr*  ■§£  +  -^]  .  (HI.6) 

In  addition  the  pressure  on  the  boundary  was  given  by 

p=p0/r3Y  (nr.  7) 

The  above  set  of  equations  were  expressed  as  finite  differences. 
The  variables  were  situated  at  different  points  in  the  finite  difference 
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lattice  to  that  the  derivatives  could  be  expressed  in  a  different  form 
accurate  to  second  order.  The  arrangement  used  is  shown  in  Fig.  9. 
The  velocity  u  was  evaluated  at  the  points  (i,  j  +  i),  the  pressure  p 
at  the  points  (i  +  i,  j).  and  the  radial  distance  r  at  the  points 
(i  +  j  +  i),  where  i,  j  are  the  coordinates  of  a  typical  point  of  the 
lattice.  At  the  boundary  the  velocity  increment  was  found  using  a  three 
point  difference  formula  for  the  derivative  8p/8y.  An  iteration  scheme 
using  averaged  values  was  employed,  in  a  manner  similar  to  that  des¬ 
cribed  in  the  previous  appendices. 

The  initial  points  were  obtained  using  Gilmore's  theory.  Given 
a  set  of  points  describing  the  motion  of  the  bubble  wall,  solutions  are 
found  along  the  outward-going  characteristics  starting  from  each  of  the 
given  points  and  finishing  on  a  certain  line  t  =  const.  The  points  along 
this  line  then  provide  the  initial  points  for  the  calculation.  The  incre¬ 
ment  in  y  between  each  of  the  points  is  found  from  Eq.(III.  1). 

Because  of  the  purposes  of  the  calculation,  the  computations  had 
to  be  performed  over  a  fairly  extensive  range  of  y.  However,  the  need 
for  accuracy  decreases  in  moving  away  from  the  bubble  wall,  and  the 
spacing  of  the  points  can  be  reduced  correspondingly.  This  is  already 
achieved  to  some  extent  by  the  initial  points  since  their  spacing  with 
respect  to  y  varies  roughly  as  r*.  Finer  accuracy  is  required 
near  the  boundary  when  the  collapse  becomes  more  rapid.  Hence  a 
scheme  of  interpolation  was  embarked  on,  similar  to  that  described  in. 
Appendix  II  for  the  method  of  characteristics.  When  the  differences  be¬ 
tween  first  and  second  estimates  at  the  boundary  become  greater  than  a 
certain  percentage,  interpolations  are  made  between  the  first 
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twenty-five  points  closest  to  the  boundary.  The  total  number  of  points 
is  thereby  increased  with  each  interpolation.  In  order  to  prevent  the  loss 
at  each  step  in  time  of  the  point  farthest  from  the  bubble  wall,  an  extrap¬ 
olation  is  made  assuming  the  liquid  to  be  incompressible.  This  assumption 
is  justified  because  the  magnitude  of  the  velocity  and  pressure  in  this  . 
region  is  small.  The  increment  in  time  was  determined  from  the  well- 
known  Courant-Levy  condition  on  the  relation  between  sound  speed  and 
mesh  speed. 

The  solutions  were  carried  up  to  the  final  collapse  point  and 
beyond  into  the  region  of  rebound,  where  a  compression  wave  forms  and 
moves  outwards,  eventually  turning  into  a  shock  front.  Beyond  this  stage 
the  solution  could  presumably  have  been  carried  on  using  methods  such 
as  those  of  von  Neumann  and  Richtmeyer^4^  or  Lax^^^  where  an 
artificial  viscosity  term  is  introduced  into  the  equations  of  motion  to 
create  a  smooth  transition  in  flow  through  the  shock.  By  this  means  the 
motion  of  the  shock  could  have  been  followed  as  it  moved  out  into  the 
liquid.  This  was  not  attempted  here,  however,  because  it  was  felt  that 
enough  information  had  already  been  obtained  to  make  am  order  of 
magnitude  estimate  of  shock  strengths  at  a  distance  from  the  center  of 
collapse;  to  have  followed  the  solution  to  any  extent  in  this  way  would 
have  involved  an  elaborate  reconstruction  of  the  numerical  procedures 
together  with  an  extensive  series  of  test  calculations. 

Up  to  the  final  collapse  point  the  Lagrangian  solution  could  be 
compared  to  that  obtained  from  the  method  of  characteristics.  Good 
agreement  was  found  once  both  programs  had  been  developed  and  tested. 
For  example  in  the  results  presented  here  the  final  collapse  radii  agree 
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to  within  lots  than  1%.  Other  comparisons  showed  a  similar  order  of 
agreement.  After  the  final  collapse  point,  the  Lag  rangian  solution 
showed  good  agreement  with  physical  considerations.  The  crest  of  the 
compression  wave  moved  with  the  local  speed  of  sound,  and  the  attenua 
tion  under  conditions  where  dispersive  effects  were  not  strong  Varied 
approximately  in  inverse  proportion  to  the  distance  from  the  center  of 
collapse. 

The  calculation  time  for  each  of  the  examples  given  here  was 


about  twenty  minutes. 
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Fig.  1  The  bubble  wall  Mach  number  ae  a  function  of  the 

bubble  radius  for  decreasing  gas  content.  The  gas 

content  is  determined  by  its  initial  pressure  pQ 

in  atmospheres.  The  index  y  has  the  value  1 .4 

and  the  ambient  pressure  p  is  one  atmosphere. 
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Fig.  2  The  bubble  wall  Mach  number  as  a  function  of  the 
bubble  radius  for  decreasing  gas  content.  The 
index  y  has  the  value  1.4  and  the  ambient  pressure 
p  is  10  atmospheres. 
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Fig.  3  The  bubble  wall  Mach  number  as  a  function  of  the 
bubble  radius  for  decreasing  gas  content.  The 
index  y  has  the  value  1 . 0  and  the  ambient  pressure 
p  is  1  atmosphere. 


Fig.  4  Variation  of  Mach  number  with  diatance  from  the 

bubble  wall  for  different  inetanta  in  time  during  the 
collapae  and  rebound  of  the  bubble.  The  initial 
internal  pre a aure  p  in  the  gaa  ia  10"9  atmoapherea. 
The  ambient  preaaure  ia  1  atmoaphere. 
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Fig.  5(a)  Variation  of  pressure  with  distance  from  the 
bubble  wall  for  different  instants  in  time  dur¬ 
ing  the  collapse  of  the  bubble.  The  conditions 
correspond  to  those  of  Fig.  4. 


Fig.  5(b)  Variation  of  pressure  with  distance  from  the 

bubble  wall  for  different  instants  in  time  during 
the  rebound  of  the  bubble. 
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Fig.  6  Variation  of  Mach  number  with  distance  from  the 
bubble  wall  for  different  instants  in  time  during 
the  collapse  and  rebound  of  the  .bubble .  The  initial 
pressure  p  in  the  gas  is  10"4  atmospheres.  The 
ambient  pressure  is  1  atmosphere. 


Fig.  7(a)  Variation  of  pressure  with  distance  from  the 
bubble  wall  for  different  instants  in  time  dur¬ 
ing  the  collapse  of  the  bubble.  The  conditions 
correspond  to  those  of  Fig.  6. 


Fig.  7(b)  Variation  of  pressure  with  distance  from  the 
bubble  wall  for  different  instants  in  time  dur¬ 
ing  the  rebound  of  the  bubble. 
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Fig.  8  Finite  difference  mesh  in  R  -  t  plane  for  solution 
by  method  of  characteristics. 
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Fig.  9  Finite  difference  mesh  in  R  -  t  plane  for  solution  of 
Lagrangian  equations. 
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